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ABSTRACT 


Current-sheet formation and magnetic reconnection are believed to be the basic 
physical processes responsible for much of the activity observed in astrophysical 
plasmas, such as the Sun’s corona. We investigate these processes for a magnetic 
configuration consisting of a uniform background field and an embedded line dipole, a 
topology that is expected to be ubiquitous in the corona. This magnetic system is driven 
by a unifonn horizontal flow applied at the line-tied photosphere. Although both the 
initial field and the driver are translationally symmetric, the resulting evolution is 
calculated using a fully three-dimensional magnetohydrodynamic (3D MHD) simulation 
with adaptive mesh refinement that resolves the current sheet and reconnection dynamics 
in detail. The advantage of our approach is that it allows us to apply directly the vast 
body of knowledge gained from the many studies of 2D reconnection to the fully 3D 
case. We find that a current sheet forms in close analogy to the classic Syrovatskii 2D 
mechanism, but the resulting evolution is different than expected. The current sheet is 
globally stable, showing no evidence for a disruption or a secondary instability even for 
aspect ratios as high as 80:1. The global evolution generally follows the standard Sweet- 
Parker 2D reconnection model except for an accelerated reconnection rate at a very thin 
current sheet, due to the tearing instability and the formation of magnetic islands. An 
interesting conclusion is that despite the formation of fully 3D structures at small scales, 
the system remains close to 2D at global scales. We discuss the implications of our 
results for observations of the solar corona. 
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1. INTRODUCTION 


Magnetic reconnection has long been postulated to be the fundamental process 
underlying most solar activity (Parker 1972, 1988; Priest 1981; Priest & Forbes, 2000). 
The basic scenario is that free energy in the form of electric currents is generated in the 
corona as a result of either the emergence of non-potential flux from below or the 
stressing of pre-existing coronal field by photospheric motions. Given the high Lundquist 
numbers (S > 10 10 ) for typical coronal parameters, these currents could persist almost 
indefinitely, but some process (or processes) brings the scale of the currents down to 
values where the frozen-in-flux condition can be broken. The magnetic free energy then 
is released to the coronal plasma by reconnection. Depending upon the properties of the 
reconnection, the energy release may result primarily in heating, as in models for coronal 
heating (e.g., Klimchuk 2006; Rappazzo et al. 2008), or in violent mass motions, as in 
models for chromospheric explosions, coronal jets, and coronal mass ejections (Lynch et 
al. 2008, 2009; Pariat et al. 2009), or even in bursts of particle acceleration, as in flare 
models (e.g.. Miller et al. 1997; Drake et al. 2006). 

It is evident from this discussion that the fonnation of current singularities, which 
we shall generically term “current sheets,” is central to reconnection and to solar activity. 
Various mechanisms have been discussed for current-sheet formation, but the one that we 
feel is most intuitive and most likely to apply to the corona is the null-point deformation 
process described in Syrovatskii’s seminal work (Syrovatskii 1971, 1978a, 1978b, 1981). 
The fundamental argument is that for a magnetic topology with discontinuous 
connectivity - in particular, at a classic 2D X-type null point (actually an X line) with 
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four topologically distinct flux systems - any stress applied to the system is likely to lead 
to discontinuous stress at the null and separatrices and, hence, to the rapid formation of 
current sheets there. This picture has been shown to be robust by many investigations, 
both analytic and numerical (Syrovatskii 1981; Antiochos 1990; Karpen et al. 1995, 

1996, 1998; Bim et al. 1998; Antiochos et al. 2002). Although Syrovatskii formulated the 
theory originally for a 2D null-point configuration, the mechanism appears to be 
physically valid in 3D as well, with only straightforward modifications (Antiochos 1996; 
Edmondson et al. 2009, 2010). 

In addition to the large-scale current-sheet fonnation process, the mechanism that 
breaks the frozen-in flux condition locally is critical for magnetic reconnection. In fact, it 
may well be that this microscale mechanism controls the energy release rate of the 
reconnection, thereby determining its global manifestations. After decades of intense 
study, a generally accepted picture has emerged for the simplest possible situation: 
steady-state reconnection in a 2D X-point topology under the MHD approximation. The 
studies have shown that if the dissipation mechanism, usually a simple isotropic 
resistivity, is spatially localized at the X point (Parker 1973; Biskamp 1993), then the 
reconnection takes on the Petschek form with a pair of slow-mode shocks, a current sheet 
short (in the 2D plane) compared to the global scale, and an inflow speed of the order of 
the Alfven speed (Petschek 1964). On the other hand, if the resistivity is roughly 
unifonn, then the reconnection adopts the Sweet-Parker (Sweet 1958; Parker 1963) form, 
with a current-sheet length of the order of the global scale and an inflow speed much less 
than the Alfven speed. The physical origin of this sensitivity to the spatial localization of 
the resistivity is well understood, discussed originally by Parker (1972) and elaborated 
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upon more recently by Kulsrud (2001). It should be emphasized, however, that both the 
Sweet-Parker and Petschek models are steady-state models only, and therefore do not 
address the problem of self-consistent current-sheet fonnation and disruption. 

In previous MHD studies of solar activity in a 2.5D null-point coronal topology 
driven by photospheric motions, it was found that the resulting evolution is well 
represented by a hybrid of the Syrovatskii and Sweet-Parker models (Karpen et al. 1995, 
1996, 1998). A smooth horizontal photospheric flow produced discontinuous stress at the 
null, resulting in its deformation to a current sheet and, hence, to reconnection there. 
Those studies all used numerical resistivity, which has no pre-defined spatial localization 
along the current sheet and, consequently, tends to produce a long global-scale current 
sheet and reconnection resembling the Sweet-Parker model rather than Petschek. A major 
new feature of these driven models, however, is that if the current sheet becomes 
sufficiently long, it is subject to tearing-mode instabilities (Furth et al. 1963; Bateman 
1978), which result in the fonnation of magnetic islands and can enhance the inflow 
speed and reconnection rate (Karpen et al. 1995, 1996, 1998; Antiochos et al. 2002). The 
general conclusion, therefore, from the 2.5D MHD studies with numerical or uniform 
resistivity is that reconnection in the solar corona is well described by the standard 
Syrovatskii-Sweet-Parker scenario, but with the addition of magnetic island formation 
and the resulting dynamics. 

A key question is whether this conclusion remains valid for a realistic 3D system 
(see Antiochos & Dahlburg, 1997). Due to the intrinsic complexity of 3D topologies, 
current-sheet formation and reconnection in such systems have begun to be explored only 
recently (e.g., Lynch et al. 2008, 2009; Pariat et al. 2009; Edmondson et al. 2009, 2010). 
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For systems that are fully 3D, such as an isolated X-point topology with no special 
symmetry (e.g., Pariat et al. 2009), the geometry is so different from the 2D model that it 
is difficult to relate the knowledge gained from the lower-dimensionality studies to the 
3D evolution. One approach to bridging the gap between 2D and 3D is to consider a 
system that is initially 2D, but to allow for a fully 3D evolution. The goal here is to be 
able to use some of our physical intuition gained from the large body of 2D reconnection 
studies in order to understand the more complex 3D system. 

Such an approach has been used to investigate the 3D evolution of a pre-existing 
ID current sheet (Dahlburg et al. 2003, 2005, 2006). Those studies produced an 
interesting and potentially important new result. As expected, the sheet tears to fonn 2D 
magnetic islands, but these are susceptible to an ideal 3D secondary instability, which 
grows much faster than the 2D tearing mode. The physical interpretation is that a 
magnetic island in 3D is actually a twisted flux rope, which can be susceptible to kink- 
type ideal instabilities that induce a burst of energy release. The results provide a physical 
mechanism for producing the type of reconnection bursts required by Parker’s nano flare 
model (Parker 1972) for coronal heating (Klimchuk 2006). Although they are interesting 
and potentially highly important, those results are not fully self-consistent, because the 
current sheet is assumed a priori rather than being generated by photospheric motions. 
Furthermore, the initial current sheet is ID and of infinite extent (in principle) in the 
transverse directions, rather than having physical dimensions determined by the global 
scales in the problem. 

In this paper we perfonn an investigation of the self-consistent formation and 
reconnection of a 3D current sheet driven by photospheric motions. The goal is to use as 
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much as possible the insights gained from the 2D studies to understand the 3D evolution. 
For the initial topology we assume a 2D X line similar to previous studies (Karpen et al. 
1995, 1996, 1998; Antiochos et al. 2002), except that we use Parker’s (1972) ansatz of 
considering the magnetic field in a Cartesian box with the bottom and top of the box 
representing photospheric regions of different polarity. As will be evident below, this 
assumption allows us to stress the field with a simple ID flow at the photosphere, which 
preserves the photospheric flux distribution there. Although the initial field is 2D and the 
driving flow is ID, the calculated evolution is fully 3D from the outset, allowing us to 
investigate the self-consistent fonnation and disruption of a 3D current sheet. The details 
of our numerical model and results are described in §§ 2 and 3 below, and we draw 
conclusions about their implications for understanding coronal evolution in § 4. 


2. NUMERICAL MODEL 


The magnetic field for our study is initially translationally symmetric and 
potential, i.e. current-free. We begin with an infinite line dipole of magnetic charge 
embedded in a uniform background field. The x component of the vector potential for this 
field is given by 


Ao ( X )= B 0 


y 


2 ay 

2 2 
y + z 


The two terms on the right-hand side represent respectively a uniform background field 
of strength Bo oriented parallel to the positive z direction, and the field due to a single line 
dipole of strength a Bo (where a > 0 is a dimensionless scale factor) positioned at ( y,z ) = 
(0,0) and oriented anti-parallel to the background field. The x direction is oriented along 
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the line dipole, and the y direction is the transverse horizontal coordinate. The general 
topology is that of a four-flux system with an X-type null point in the yz plane (see Fig. 

1). The separate flux systems are comprised of field lines that close from within the line 
dipole to the lower z boundary to either side in y or, alternatively, that connect the upper 
and lower z boundaries to each other. 

Figure 1 also indicates schematically the footpoint motion that is imposed at the 
top of our Cartesian box. To allow for an arbitrary horizontal displacement of the 
footpoints at the top (z) boundary while simultaneously avoiding the complexity of 
inflow and outflow conditions at the side (y) boundaries, we augmented the magnetic 
field above to make it periodic along the y direction. An infinite series (N—*o o) of pairs of 
additional line dipoles positioned at z = 0 andy = ±nd for n . = 1,2, . . ., N, where d is a 
prescribed periodicity length, contributes the following to the vector potential: 



Although the sum is strictly periodic only as /V— >oo, we find in practice that N = 25 terms 
are numerically sufficient on our domain. The resultant initial magnetic configuration is 
shown in Figure 1 over the periodic y interval [—d/2 ,+d/ 2] and the vertical z interval 
[h,h+d/ 4], where we have chosen h = d/ 40. For completeness, the convergence of the 
infinite series is demonstrated in the Appendix. 

There remains the specification of the x extent of our domain. Again, to minimize 
any influence of the boundaries, we elected to impose periodic conditions along the 
initially translationally symmetric x coordinate. Anticipating that y displacements of the 
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order of the height of the domain, Az = d/4, would be imposed and introduce structure on 
that or smaller scales, we chose the x domain extent to be |jc| < d/4. 

We calculated the evolution of the field in Figure 1 with the Adaptively Refined 
Magnetohydrodynamics Solver (ARMS) (Welsch et al. 2005; DeVore & Antiochos 2008; 
Lynch et al. 2008, 2009; Pariat et al 2009; Edmondson et al. 2009, 2010). We solved the 
following set of 3D compressible, ideal MHD equations: 

^ + V*pv= 0 
dt 

d -^ + V • pvv + VP = — fvxfi 
dt 4 Jt ' ' 

'I rj-1 

+ vrv + (y-2)rv«v = 0 

dt 

— = Vx(vxB ) 
dt ' ' 

where all variables have their usual meanings. In the adiabatic equation for the 
temperature, the ratio of specific heats is y = 5/3. We have used the ideal gas law P = pRT 
for the plasma equation of state, where R is the gas constant. Gravity is not included in 
this calculation, since we are interested in the local structure and evolution of current 
sheets in the high-temperature corona. 

The MHD equations above can be made dimensionless by extracting scale factors 
for the mass density p s , magnetic field strength B s , and characteristic length L s . The 
characteristic length L s , for example, can be set to that of a coronal structure of interest; 
the associated normalized length in our simulation is h, the depth of the line dipole. The 
mass density p s and field strength B s are similarly free parameters; their nonnalized 
counterparts po and B 0 will be specified below. As a result of this scaling, all values of 
the MHD variables in the simulation are of order unity. For example, the initially 
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uniform mass density po and background field strength Bo are chosen to set the global 
average Alfven speed V A (defined below) to unity. We set the relative line-dipole strength 
a so that the null is positioned well inside the simulation domain, insulating the current- 
sheet formation and dynamics from any uncontrolled boundary influences. The initially 
uniform thennal pressure po is chosen to assure that the low-beta condition (thennal 
pressure much smaller than magnetic pressure, so that the dynamics are magnetically 
dominated) holds everywhere in the domain other than in the immediate vicinity of the 
null line. Figure 2 illustrates the plasma beta across the inner portion of the magnetic 
structure, the white contour showing the [3 = 1 surface. Only within this small region in 
the neighborhood of the null line does the plasma beta rise substantially above unity. We 
emphasize, however, that although the system as a whole is low-beta, the thermal 
pressure nevertheless plays an important role in the evolution. It dominates in the small 
region near the null line, so that the dynamics of current-sheet formation and subsequent 
reconnection are critically dependent upon the thennal pressure being finite rather than 
zero. 

In the low-beta limit, driving motions that are slow compared to the Alfven speed 
yield an approximately force-free evolution. This regime is appropriate for current-sheet 
fonnation in the corona in response to slow photospheric driving. Much like the plasma 
beta, the Alfven speed varies most strongly over a small region in the vicinity of the null 
line, where it approaches zero. Elsewhere, the Alfven speed varies smoothly, approaching 
a uniform value well away from the line dipole and its periodic images. It is convenient 
to define the global average Alfven speed (V A ) = ^ 2 E u / M , where E M is the total 
magnetic energy and M the total mass in the system. The scale factors for mass density 
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(p s ) and magnetic field (B s ) can be chosen such that <Vj> = 1 in our initial setup in 
Figure 1. As a consequence, all velocities in our simulation will be measured in units of 
the initial global average Alfven speed < V A > of the corresponding coronal configuration. 
We noted earlier that the length scale factor L s fixes the scaled line-dipole depth h, which 
defines a global average Alfven time (t A ) = h/(V A ). We can set h = 1 , whence <t A > = 1 . 

As a result, the elapsed time in our simulation similarly will be measured in units of the 
initial global average Alfven time <t A > in the corona. 

We energize the system by imposing a sub-Alfvenic, spatially uniform flow all 
across the upper z boundary, as indicated schematically in Figure 1. In order to minimize 
transient effects as the motions get underway, we employ a gradually accelerating cosine 
profile in time: 

V-V t (l)y 




K 


I 


1 - cos 


71 — 


\ 


0 <t<t n 


[ K t>{ 0 

The magnitude of the imposed flow is Vo « 1, slow compared to the global average 
Alfven speed, while the duration of the smooth acceleration phase is to » 1, long 
compared to the global average Alfven time. An accumulated displacement IVo/2 ~ h 
distorts the null line sufficiently to create a substantial current sheet by the end of the 
acceleration phase. After time t = to, the driver speed is held constant, so that any 
subsequent burstiness in the system will be an intrinsic property of the dynamics at the 
sheet rather than a response to a time-variable driver. Line-tied conditions are imposed at 
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both z coordinate boundaries, with the footpoints held fixed (zero flow speed) at the 
lower boundary. Consequently, the only way that the magnetic field can change topology 
is through reconnection. 

Figure 3 shows the initial numerical grid for the simulation. Our base grid 
consisted of 2x2x1 blocks distributed uniformly in {x, y, z}, with 8 3 grid points per block 
(note that only every other grid point is shown in the figure). We then imposed an initial 
refinement to a maximum of 6 levels, over a subvolume that we anticipated would be 
sufficient to encompass fully the developing current sheet. This region, which has the 
highest density of grid lines {cyan) in Figure 3, is centered rightward and upward from 
the null line. Beyond this subvolume, the grid coarsens by a factor of two at each change 
in refinement. We employed adaptive mesh refinement (AMR) to increase further the 
grid resolution near the strengthening current sheet as the simulation progressed. Our 
mesh-refinement criterion tests the gradient of the magnetic-field magnitude in regions of 
weak to moderate field strength. The grid was allowed to refine up to a maximum level of 
9, corresponding to a finest unifonn grid of 4096 x 4096 x 2048 cells if employed over the 
entire domain. The resultant minimum grid spacing, 8 = 2.44x 10’ 4 d, is far smaller than 
the scale of the initial magnetic structures. 

3. RESULTS 

3.1. Initial Configuration and Current Sheet Development 
The domain size is set to d = 40/z = 40, so that the simulation was performed on the 
space {(. x,y,z ) | x G [-10,10] ,y G [-20,20] , z G [1,11]}. We chose a relative line-dipole 
strength a = 20/3, which places the initial null line near domain center in the z direction 
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(see Figure 1). The scale factor B s for the absolute field strength is chosen so that the 
simulation background field Bq = 3. We found that the root mean square field magnitude 

in our initial configuration is approximately Vt# . Choosing po= 1 then yields the global 
average Alfven speed <V A > = 1, as desired. We set the thermal pressure pn = 0.01. The 
global average plasma beta therefore is <f3> = 0.02, while the minimum value (3 min ~ 10~ 4 
is attained at the bottom 2 boundary, in the closed field of the line dipole (see Figure 2). 

For the imposed flows at the top 2 boundary, we chose a peak speed Vo = 0.045 and 
a corresponding duration for the acceleration phase to = 100. The resultant time profile of 
the boundary flow speed is shown in Figure 4. At times t>to= 100, the flow is held 
steady at its peak value. Since the plasma and field are fixed at the bottom 2 boundary, we 
can expect the average velocity over the entire volume to approach roughly Vo/2. For the 
assumed values of po, Vo, and system size d, the predicted asymptotic kinetic energy 
integral over the volume is Ek = 2. The simulation closely adheres to this expectation, as 
shown by the time profile of kinetic energy in the top panel of Figure 5. 

In the early development stage, t < 40. 19, the driving speed and footpoint 
displacements ramp up slowly and only weak currents appear in the corona. If the line 
dipole were not present, so that the initial field was uniform and vertical everywhere, then 
the field would remain nearly potential throughout the evolution. The photospheric 
motions would generate a uniform B v that increases linearly with the total photospheric 
displacement in the y direction, and the field lines simply would tilt increasingly away 
from the vertical as their top footpoints suffer ever greater displacements. The closed 
bipole flux, however, acts like an obstacle to the background field, and the field must 
deform around the bipole flux as it tilts. The field strength is small in the neighborhood of 
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the null, so that field lines passing through that region cannot sustain large stresses. Those 
above the null connect to the top z boundary, and the imposed flows displace them 
laterally with relative ease. In contrast, those below the null connect to the bottom z 
boundary, where the field magnitude is quite large and the footpoints are fixed. Inside the 
bipole near the line source, strong volumetric currents appear as the field resists the 
stresses transmitted by the motions imposed at the top boundary. 

The process of current-sheet formation is intuitively simple. Prior to onset of the 
boundary driving motions, inner and outer spine lines (actually planes in 3D) depart the 
null point toward the bottom and top z boundaries, respectively, and are perfectly vertical. 
The imposed driving motion moves the top footpoint of the outer spine to the right, which 
imparts a stress to the field and drags the lower end of the outer spine to the right. The 
inner spine, however, feels no stress and stays more or less fixed. Consequently, the 
spines dislocate, causing the null point to defonn. This process begins as soon as the 
motions get underway, initiating the development of the current sheet. The small but 
finite gas pressure limits the narrowness of the current sheet at the outset. Figure 6 shows 
the system configuration at the initial and final times, t = 0.0 and 147.64 seconds, (the 
corresponding movie is available in the electronic version of the paper). The field lines 
are shown in white against the perpendicular current contour background (the strong 
current, \J X \ > 5, region is cyan). The singular current sheet is clearly seen between the 
dislocated spine lines in the left figure. 

As the driving flows accelerate, the null line stresses further, fonning an 
increasingly strong neutral sheet. The adaptive mesh refinement first turns on at time t = 
40.19 (see Figure 5), as the null defonns substantially and its current structure steepens to 
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approach the grid scale at the initial maximum refinement level (6). At the successive 
times t = 40.19, 50.317, and 60.187, the maximum refinement level increases to 7, 8, and 
finally 9, as the grid subdivides to resolve the continually thinning sheet. The ongoing 
boundary motions further lengthen the current sheet and increase the volume occupied by 
grids at the highest refinement level (9), as illustrated at time t = 70.003. At each of these 
discrete times, as well as at other intervening times in the simulation, there are strong 
changes in the slope of the kinetic energy plot. These evidently are associated with the 
internal restructuring of the current sheet as the grid refines in discrete jumps. Through 
the remainder of the acceleration phase and during the subsequent phase of constant 
imposed boundary motions, at times t > 70.003, the grid has reached its maximum 
refinement level (9) over the entire sheet between the inner and outer spine lines. At this 
level of refinement, the current sheet half-thickness a is approximately the finest grid 
spacing, 2.5><1 O' 4 d= 1.0x1 O' 2 . 

At times t <70, the current sheet lengthens and intensifies as the driving flows 
accelerate. Subsequently, the lengthening subsides as the imposed flows reach top speed. 
At times t > 100, the driving flows are constant. During this late phase, the current sheet 
exhibits a more or less steady state of approximately fixed length (see Figure 7), in which 
plasmoids repeatedly fonn within the sheet and are ejected out of one end or the other, as 
described in detail below. In order to interpret our findings using results from MHD 
reconnection theory, it is necessary to estimate the effective magnetic diffusivity rj in the 
simulation. Balancing the inflow speed V, of magnetic flux into the current sheet against 
the resistive annihilation across the sheet yields 

V L= rj_ 

2 ’ 

a a 
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whence 


17 = V t a . 

Combining this estimate with a characteristic global scale L, such as the half-length of the 
current sheet, and the local Alfven speed V a, the effective Lundquist number S then is 

rj V ]a 

Averaged along the current sheet for times t > 100, we find that the current-sheet inflow 
and Alfven speeds are approximately V, = 0.03 and V a = 0.45, respectively. The spine-to- 
spine separation length is 2 L = 0.04J = 1.6, and as noted above the half-thickness of the 
sheet is a = 0.01. Thus, the effective numerical resistivity and Lundquist number for the 
fully developed, high-resolution current sheet can be estimated to be rj = 3x 10 " 4 and S = 

1 ,2x I 0 3 , respectively, during the approximately steady-state phase of the evolution. It 
should be emphasized, however, that these values for rj and S apply only at the current 
sheet, itself. Outside the current sheet, the magnetic structure is very well resolved and 
the evolution is effectively ideal. 

3.2 Plasmoid Formation and Ejection 
The relatively long-wavelength resistive tearing mode is the most important MHD 
instability in the developing current sheet. In 2D, resistive tearing of the sheet generates 
alternating X- and O-type field topologies all along its length (Biskamp 1993). Flux 
accumulating around the O points creates magnetic islands of locally closed field lines, 
which trap the reconnection-jet plasma and thus form enhanced-density plasmoids. In 3D, 
the additional degree of freedom enables the development of much richer topologies. In 
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the absence of a guide field, as in our case, 3D perturbations with variations along the 
nominal symmetry axis generate quasi-ergodic, slinky-like field lines via reconnections 
between neighboring cross-sectional planes. Current sheets of finite length, whether 2D 
or 3D, eject their fully developed magnetic islands/plasmoids in the direction of the 
closest tenninus of the sheet. Therefore, we differentiate the tearing-mode process into 
two stages: plasmoid fonnation and plasmoid ejection. The formation stage consists of 
the initial field line pinching and linear growth of the plasmoid; its characteristic length 
scale is the current sheet half-thickness (i.e., a = 0.01). The ejection stage occurs once the 
plasmoid has grown to the order of the current-sheet thickness, so that nonlinear effects 
come to dominate the dynamics. We discuss the final scale of the ejected plasmoids in the 
following section. 

The linear growth time r of the tearing mode is intermediate between the short 
Alfven time r A = a/V A and the long resistive time r, = a 2 /?/ of the sheet (Furth et al. 1 963). 
It depends upon the large-scale properties of the current-sheet magnetic field, the sheet 
thickness, the wavelength along the sheet, and a weighted geometric average of the two 
characteristic time scales. The linear growth rate can be expressed as (e.g. Bateman 1978) 


Y ~ 0.9 


(kaf 5 (A 'of 5 



where A’ is the jump in the logarithmic gradient of the magnetic-field perturbation across 
the current sheet and k is the wavenumber along the sheet. The quantity A’ depends solely 
on the properties of the ideally evolving field, and at moderately long wavelengths ka < 1 
can be estimated as (Furth et al. 1963) 
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The tearing-mode growth rate then becomes 


1.6 


( kaY * 


t/ 5 t /5 

A rj 


which can be inverted to give the growth time 

r » 0.6 (kaY 5 

previously used to interpret tearing and island fonnation in a strictly 2D current sheet 
(e.g. Karpen et al. 1995). 

Throughout the evolution of our 3D system, the discernible tearing wavelengths 
are between one-fourth and one-half of the total length. For a current-sheet aspect ratio 
L/a = 80 and at most two wavelengths along the sheet half-length L, we find 

2 / 3 / 

r < 0.3r ./X, • 

Recalling that the current-sheet Alfven speed V a = 0.45 and effective magnetic diffusivity 
rj = 3x 10" 4 , we find r A = 0.02 and r„ = 0.3. The linear growth time during the plasmoid 
fonnation stage then is r< 0.03. This is much shorter than the characteristic inflow time 
for the sheet, r, = a/Vi = T, = 0.3, so that applying the linear perturbation analysis for an 
equilibrium sheet to our simulation is reasonable. In addition, many e-folding times 
elapse during the late sheet-lengthening and steady-state phases of the evolution. This 
provides ample opportunities for new islands to be seeded at small amplitude, grow to 
discernible size, and saturate in the nonlinear regime. For example, the top panels of 
Figure 8 (at times t = 1 1 1.66, 1 12.02, and 1 12.54) illustrate emergence from the noise of 
a sheet-spanning plasmoid over some thirty linear growth times. 

For the most part, we find that the plasmoids initially form near the center of the 
current sheet. Since the typical tearing wavelength is one-fourth to one-half of the sheet 
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length, we find on average that there are two to four irregularly shaped plasmoids in any 
given perpendicular vz plane. The plasmoids generally are well separated from one 
another along the sheet within the vz plane. Nonlinear effects saturate the rapid growth of 
the plasmoid when its width becomes comparable to that of the current sheet: at that 
point, the locally available free energy of the magnetic field is exhausted and a new 
quasi-equilibrium state is attained. Thereafter, the plasmoid evolution is governed by the 
gradually evolving state of the current sheet as a whole. 

The ejection time for individual plasmoids is somewhat variable during the 
evolution, depending upon the precise location of formation along the current sheet and 
the instantaneous length of the sheet itself. The ejection speed reflects the local outflow 
velocity of the reconnection jets within the sheet, which is approximately the Alfven 
speed V a = 0.45. Plasmoids sweeping out the ends of the sheet - especially the larger 
plasmoids - split the sheet there, temporarily reducing its length. A characteristic time 
scale x e for plasmoid ejection is the Alfven time across the maximum sheet half-length L : 
r e = L/Va = 1.8. Figure 8 illustrates a typical plasmoid ejection in the panels at times t 
between 1 13.08 and 1 14.95, during which a plasmoid that fonned roughly in the center of 
the sheet is ejected from it over the prescribed Alfvenic ejection time scale r e (the 
corresponding movie is available in the electronic version of the paper). 

Our simulation is fully 3D at all times, and there is no guide field imposed along the 
initially symmetric axis of the magnetic configuration. Consequently, the tearing-mode 
perturbations that seed the formation of the magnetic islands and plasmoids in each yz 
plane do not necessarily develop coherently over any significant length along x between 
the planes. In this simulation of a truly neutral sheet, the only guide-field component that 
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emerges is seeded by 3D numerical perturbations. This limits the correlation length along 
the nominal symmetry axis of the current sheet to only a few grid cells. 

As a consequence of this lack of correlation along the third dimension, together 
with the short time for island growth relative to the long time for sheet development, the 
plasmoids that form are short and evolve effectively independently on different yz planes. 
Although there is some evidence for kinking of these short structures, they do not become 
sufficiently large to affect the overall integrity of the sheet. We illustrate this point in 
Figure 8 (and the corresponding movie available in the electronic version of the paper), 
which shows the current sheet in the xy plane as viewed from above. The times are the 
same as those in Figure 8, and the white stripes mark the vz plane shown in cross-section 
in the previous figure. Purple shading shows strong current-density magnitude (M> 5), 
while red shading highlights the plasmoid density enhancements (p> 2). The white 
arrows in Figure 8 track the plasmoid formation, evolution, and ejection indicated by the 
black arrows in Figure 8. The 3D plasmoids evolve quasi-independently along the third 
dimension, at some times separating and at other times coalescing as they interact during 
their transit along the sheet. Although the resultant structures occasionally kink, they 
never extend very far along the third dimension and are ejected quickly enough that the 
current sheet always maintains its overall stability. It is obvious from Figure 8 that the 
plasmoid structure is highly irregular and fully 3D at all times. 

3.3 Quasi-Steady Reconnection Dynamics 

The late-phase, approximately steady state of our system ( t > 100) consists of a 
dynamically stable, 3D neutral current sheet with an aspect ratio of around 80:1. Clearly, 
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the system dynamics favor a rather quiet steady-state reconnection with occasional 
fonnation and ejection of localized magnetic islands, rather than an explosive or bursty 
reconnection that results in large-scale turbulent breakup of the current sheet. These 
dynamics are a direct consequence of the nature of the driving flows, which produce 
magnetic stresses localized primarily to the separatrices between flux systems. These 
stresses relax directly via the observed reconnection. To that end, the fully developed 
current sheet attains a more or less constant length between the displaced spine lines, 
although the instantaneous length in any particular yz plane fluctuates somewhat with 
each plasmoid ejection. This implies that the field-stressing and reconnection rates have 
converged to a common value, at least in a time-averaged sense. 

Figure 10 shows the color-shaded velocity magnitude at the x = 0 mid-plane for 
several times t = 100.1, 1 1 1.66, 1 15.74, and 147.64, typical of the phase of statistically 
steady dynamics. The high-speed (\V\ > 0.2) reconnection jets are clearly seen as the 
purple-shaded regions within and at the ends of the current sheet. White arrows show the 
vector velocity direction, essentially perpendicular to the sheet at either side of it, and 
principally along the sheet within it. Black contours illustrate the extent of the strong 
currents (|J| > 5) in the system. Both the average inflow speed V, across the sheet and 
average outflow jet speed V 0 along the sheet increase gradually with time. In addition, the 
flow speeds below the current sheet are slightly lower while the Alfven speeds are 
slightly higher, compared to those above the sheet. This reflects the higher field strength 
within the closed line-bipole region, which is more resistant to stressing. Table 1 
enumerates the average inflow and outflow speeds and their ratios, at the times shown in 
Figure 10. We note that there also is a slight rising trend versus time of the reconnection 
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rate, as measured by the inflow/outflow ratio. Evidently, the sheet has not yet attained a 
perfectly balanced steady state, as it continues to adjust to the fixed driving motions and 
the intennittent dynamics of plasmoid formation and ejection. 


It is instructive to compare our results with those predicted for steady-state Sweet- 
Parker and Petschek reconnection (e.g. Kulsrud 2005). The values in Table 1 show that 
the average ratio of inflow to outflow speeds over the interval is VJ V 0 = 0.11, whereas the 
ratio of current-sheet thickness to length is a/L = 0.012. These numbers are very different 
from each other and from the common value predicted by Sweet and Parker, 1 Hs = 0.03, 
for our Lundquist number S = 1 ,2x 1 0 3 . Compared to Sweet-Parker, our quasi-steady 
current sheet is thinner yet its inflows are faster, by factors of about three. In steady-state 
Petschek theory, on the other hand, the reconnection region along the current sheet 
shrinks to an effective half-length L' given by 
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For our Lundquist number S, the predicted inflow/outflow ratio VJ V 0 = 0.19, and the 
reconnection region shortens to L’= 0.02 so that a/L' = 0.5. The latter is commensurate 
with the size of our magnetic islands; however, the Petschek inflow speed and associated 
reconnection rate are too high by nearly a factor of two. Thus, our dynamically tearing 
current sheet exhibits an overall Sweet-Parker-like structure with locally reconnecting 
regions at the Petschek scale, but with an inflow/outflow ratio and reconnection rate that 
are intennediate between these 2D, steady-state reconnection models. 
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The key point is that our driven current sheet is physically different than that of 
the Sweet-Parker or Petschek reconnection models. In these steady-state models the 
reconnection rate, (the inflow velocity), is the important derived parameter; but, for our 
current sheet, the inflow/outflow ratio is fixed by the initial and boundary conditions. The 
current sheet forms due to the horizontal displacement (at imposed speed Vo) of the 
initially uniform vertical field (at strength B o and mass density po). Thus, the inflow 
speed into the sheet is V t ~ Vo = 0.045 and the outflow speed within the sheet is the local 
Alfven speed, V 0 = Va ~ B 0 /2^4npo = 0.42, using Bo/2 as an average field strength within 
the null region. Their ratio VJV 0 = 0.1 1 is fully prescribed by these parameter choices; 
hence, the question of the reconnection rate is not relevant to our system. 

The important question that does remain, however, is the scaling of our current 
sheet scale with increasing S, which is equivalent to increasing grid refinement. There are 
four critical parameters that describe a general reconnecting current sheet. The first and 
most important is the current sheet width, a, or more exactly, the scale at which the 
frozen-in flux condition is broken and the reconnection of field lines occurs. For our 
calculations, which use numerical resistivity, the flux-breaking scale is simply the grid 
size. Although this choice may seem arbitrary, it is actually the most appropriate 
resistivity model for an MHD calculation. The flux breaking in the Sun’s corona and in 
most astrophysical plasmas is expected to occur at kinetic scales such as the ion skin 
depth or gyroradius; therefore, no MHD model can capture the physics of the 
reconnection region rigorously. All the kinetic studies have shown, however, that the 
reconnection sets in at some fixed scale, such as the ion skin depth or gyroradius, (e.g., 
Cassak et al 2005; Yamada 2007). Therefore, the best that can be done in an MHD 
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calculation is to set the diffusion to occur at as small a fixed scale as possible. Our ARMS 
code does exactly this by adding diffusion to the ideal MHD solution whenever the 
current structure gets down to the grid scale. 

The next key parameter in the problem is the length of the current sheet, L. In our 
system, this length is determined by the current sheet formation process and the scales of 
the initial magnetic structure. The scale of the high-beta region around the initial null 
point is of order the global scale. In order for this region to deform to a neutral sheet 
whose width decreases down to the grid scale, the region must suffer a significant 
deformation, at least, of order its initial size. This simple argument implies that the length 
of current sheets that form in the Sun’s corona must be of order the global size of the 
relevant magnetic flux systems, as found in our simulation. 

The final two parameters are the inflow and outflow velocities. The former is set 
directly by the imposed boundary flow. The latter is set by the Alfven speed inside the 
embedded bipole. Note that the current sheet fonns, because the stress added to the 
external field causes it to press down against the field of the bipole. It does not require an 
enormous defonnation of the bipolar field to decrease the current sheet width down to the 
grid scale; hence, the field in the inflow region just outside the current sheet is of order 
the initial bipole field. In other words, there is no significant flux pile-up in our system. 

We conclude that the scale and magnitude of the global fields and imposed 
driving motions completely determine all the important parameters of the reconnecting 
current sheet. But this result seems impossible, because there are certain constraints that 
the system must satisfy, such as mass and flux conservation, which should impose 
additional, independent constraints on the current sheet parameters. The existence of 


24 



these additional constraints is the reason for the appearance of magnetic islands 
(plasmoids). The plasmoids break up the global-scale current sheet into substructures that 
introduce new scales to the system; thereby, allowing the conservation constraints to be 
satisfied while the global current sheet parameters remain fixed by the input scales and 
driving. 

Consider, in particular, mass (or flux) conservation. For the simplest possible 
current sheet structure, the Sweet-Parker case, this constraint requires that VJVi =L/a; but 
we find that V 0 /Vi= 9.1; whereas, L/a = 83. The critical point, however, is that the mass 
flux out of the global current sheet layer is carried not by the individual reconnection jets 
that have width, a, but by the expulsion of plasmoids. These have an average width, w, of 
order 10, which is sufficient to balance the mass flux into the global sheet. Note that this 
result implies that the plasmoid width, w, should be independent of Lundquist number, or 
equivalently, the grid scale of the simulation. Furthermore, the plasmoid width will scale 
directly with the driving velocity, V,, assuming that the initial conditions of the simulation 
are unchanged. 

Mass conservation also determines the average number of plasmoids in the global 
current-sheet layer. As the plasmoids separate, the current sheet between them grows to 
the maximum length compatible with mass conservation, Lm = a VJVi, beyond which 
new plasmoids must appear. For fixed driving and Alfven speeds, this limit implies that 
the number of plasmoids scales with the current sheet width, a, which for our current 
sheet is proportional to the Lundquist number. On the other hand, the maximum size of 
the plasmoids is independent of Lundquist number and set only by the reconnection rate 
and the global current sheet length. We conclude, therefore, that driven current sheets in 
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the Sun’s corona will develop a distribution of plasmoids ranging in scale from formation 
at a scale of order a Vo/Vi, to expulsion at a scale of L V/V 0 . This result is qualitatively 
similar to what is found for studies of current sheets initially in equilibrium, such as a 
Harris sheet (e.g., Loureiro et al. 2007; Bhattacharjee et al. 2009). The exact scaling of 
plasmoid number with S may be somewhat different depending on whether one is 
examining a driven or an unstable system, but all the studies appear to agree that the 
density of plasmoids increases as some power of Lundquist number. 

We note finally that the fraying and splitting at the ends of the current sheet in 
response to plasmoid ejections create localized, transient Petschek-like geometries there. 
Figure 1 1 illustrates this effect by comparing the strong current magnitude and plasma 
density, together with overlying magnetic field lines. In the bottom panel, the field lines 
(i white ) clearly exhibit sharp Alfvenic bends where the current sheet splits at the right 
side. These Petschek-like structures at the ends of the sheet are very short and transient, 
developing as the current sheet splits when plasmoids are ejected and dissipating as a 
single sheet refonns via the subsequent relaxation of the field. 

4. Conclusion 

The calculations presented in this paper constitute, in a sense, a minimalist extension 
of 2D current-sheet formation and reconnection to the fully 3D regime. The advantages 
of our approach are that it lends itself to straightforward comparisons with previous 2D 
and 3D work and allows us to employ very high numerical resolution at the current sheet. 
The results have a number of important implications for solar activity and for magnetic 
reconnection more generally. 
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First, it is clear that a current sheet forms due to the imposed photospheric 
motions, even though the motions are the smoothest possible uniform flow. This result 
demonstrates convincingly that a multipolar topology, even the extremely simple one of 
our model, inevitably will lead to current sheets in the corona. The fonnation process is 
physically identical to the Syrovatskii model: an initial X-type null defonns under the 
applied stress into a neutral current sheet bounded by two Y-type nulls. For our particular 
magnetic configuration and photospheric driving, the currents are confined to the 
deformed null-point region and are small at the separatrix surfaces. If we were to apply a 
flow at the upper surface that drives the field in the initially symmetric direction, parallel 
to the x axis, however, current singularities would form along the separatrices as well. 

The current-sheet formation process is rapid and efficient. It is evident from 
Figure 5 that the structure at the deformed null thins down to the maximum resolution for 
a fairly small displacement of the outer spine, less than 20° of tilt from the vertical. We 
expect the thinning to occur exponentially (e.g. Antiochos et al. 1999), so that even for 
the very large Lundquist numbers of the real corona, current structures down to the 
dissipation scale form on essentially ideal time scales. An important implication of this 
result is that current sheets are likely to be present in the corona even though the 
magnetic field there is generally observed to be near the potential state. We see from 
Figure 5 that a modest distortion from the potential state is sufficient to generate structure 
at the dissipation scale. Note, however, that these current singularities occur at very 
specific locations: the topological boundaries between different flux systems. It remains 
an open question whether current sheets form ubiquitously throughout the corona (Parker 
1972; van Ballegooijen 1985; Antiochos 1987). 
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Once the current sheet lengthens sufficiently that reconnection can compete with 
the slow driving motions, the sheet stops growing. From this point on (Figure 7), the 
structure near the null is quasi-static, with the stretching by photospheric displacements 
balanced by the shrinkage due to reconnection. Of course, the system is not in a true 
steady state, because the flux in the right-hand arcade of Figure 7 decreases continuously 
due to the reconnection, while the flux in the left-hand arcade increases. Eventually the 
flux imbalance will reach the point at which the reconnection stops, but we see no 
evidence for this saturation, at least for the footpoint displacements in our simulation. 
Although not extreme, the displacements clearly are substantial. Figure 7 shows that by 
the end of the calculation, the footpoints on the upper boundary have moved a distance on 
the order of the size of the closed bipole. It is intriguing that in this quasi-steady state, the 
angle of the uniform background field with respect to the vertical is roughly 20°, close to 
the stress angle required by Parker for his coronal-heating model (Parker 1983). We 
emphasize, however, that the topology of our system is significantly different from that 
of Parker’s. 

A striking result of the simulation is the extreme stability of the current sheet and 
reconnection dynamics, as shown in Figures 6 and 9. The dynamics are much less 
explosive than our simulations of reconnection with isolated 3D null-point topologies 
(Pariat et al. 2009; Edmondson et al. 2009, 2010), and even less impulsive than the 2.5D 
simulations with null lines (Karpen et al. 1995, 1996, 1998). The origin of the differences 
between the results here and our previous simulations of coronal reconnection lies in the 
fonn of the driving. First, the photospheric motions in the present calculation deform 
only the large-scale external flux system, not the small-scale embedded bipole; 
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consequently, there is no need for the free energy to spread outward from an isolated 
small region. Second, the motions inject no helicity into the corona; the field lines are not 
twisted or sheared. Therefore, reconnection, by itself, can lower the magnetic energy 
toward the potential-field minimum without having to propagate helicity over global 
scales. The implications for the corona are that slow photospheric driving of the large- 
scale coronal field can result in quasi-steady energy release via reconnection, without any 
dramatic dynamics. Such a process is exactly what is required to account for the near- 
steady heating of solar quiet regions. 

The result that the current sheet is so stable is especially surprising, given that the 
sheet achieves aspect ratios as high as 80:1. The dynamics are generally consistent with 
steady-state Sweet-Parker reconnection except for the addition of magnetic-island 
fonnation. These plasmoids carry most of the mass and are ejected from one current sheet 
terminus or the other on the Alfven time scale for the sheet length. We emphasize, 
however, that structure in the x direction does fonn at both small and global scales, and 
that this 3D structure is critically important. The small-scale structure is due primarily to 
the formation of plasmoids at random locations along the initial symmetry axis of the 
sheet. The large-scale structure is evident as a periodic arrangement of the field lines with 
exactly four wavelengths across the domain (see Figure 12). Clearly, this structure is 
seeded by the finite grid of the calculation, as is to be expected, but only the small-scale 
features are physically significant. 

We claim that the small-scale structure in the x direction is directly responsible 
for the lack of burstiness in the reconnection. Figure 8 shows that with regard to the 
plasmoids, there is absolutely no symmetry in the x direction. The plasmoids are fully 3D 
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and form at random locations throughout the current sheet. This 3D structure suppresses 
the secondary instability that generated the bursts of reconnection seen in previous 
current-sheet studies (Dahlburg et al. 2003, 2005, 2006). In those simulations, the system 
was seeded with 2D perturbations that produced magnetic islands (actually twisted flux 
tubes, owing to the finite guide field) with infinite length in the x direction. The 
secondary instability is basically an ideal, kink-like mode of the highly twisted flux tubes. 
Two-dimensional perturbations of a 3D system, however, are somewhat unphysical, 
because they assume that the noise has very long coherence length along a preferred 
direction. In the present simulation, the noise is fully 3D and leads to plasmoids whose 
coherence scale length in the x direction is similar to that along the current-sheet length in 
the yz plane. Hence, there is essentially no energy available to be released by a kinking 
instability. Furthermore, the 3D nature of the magnetic islands that we find is likely to 
limit the effectiveness of the islands in trapping particles, which could have important 
implications for flare acceleration mechanisms (e.g. Drake et al. 2006). 

An important caveat to these conclusions is that there is no guide field in our 
simulation. Due to the fonn of the driving, a true neutral sheet is generated. It may well 
be that if a guide field were present it would introduce a much longer coherence scale 
length along the x direction, so that much longer flux-tube-like magnetic islands would 
form. The situation is analogous to the issue of flare -ribbon formation. A strong guide 
field provides the reconnection coherence required to produce the distinctive ribbon 
structure. On the other hand, if the guide field is too strong, it will suppress any kinking 
via secondary instability (Dahlburg et al. 2003, 2005, 2006). We conclude, therefore, that 
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simulations with varying amounts of guide field should be the next step in studying the 
fonnation and reconnection of coronal current sheets. 
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Appendix A 

To demonstrate that the vector-potential series 



converges, note that the last term in the denominator is positive definite and also that the 

* 2222 * 2 

quantity e = (y +z )/d ' satisfies 0 < e < 1 on our simulation domain. Hence, 
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which is finite. Thus, the original series s(y,z;d) is convergent. 
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Time t 

100.1 

111.66 

115.74 

147.64 

Average Inflow Speed V t 

0.024 

0.033 

0.034 

0.043 

Average Outflow Speed V a 

0.237 

0.302 

0.278 

0.349 

Ratio Vi/Vo 

0.102 

0.109 

0.122 

0.123 


Table 1. Inflow and outflow speeds, averaged along and across the current sheet, and 
their ratio, for the four times shown in Figure 10. 
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Figure 1. Initial magnetic field, consisting of a four- flux system in the yz plane 
perpendicular to the translational symmetry axis (x). A spatially uniform driving flow is 
imposed at the top boundary. 
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Figure 2. Plasma beta. Magnetic field lines are color-shaded according to the beta values 
along them: red is low, purple is high. The white contour marks the j3= 1 boundary. 
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Figure 3. The initial grid in the perpendicular (yz) plane. Grid lines {cyan) show every 
other grid cell, up to refinement level 6 (finest) encompassing the null line. The grid is 
invariant in the translationally symmetric (x) direction. 
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Figure 4. Time evolution of the driving flow speed and, for comparison, the initial global 
average Alfven speed. 
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Figure 5. Time evolution of integrated kinetic energy (top) and of current-sheet intensity 
and extent, along with grid refinement (bottom). The adaptive mesh refinement refines 
the grid from level 6 to levels 7 through 9 at times t = 40.19, 50.317, and 60.187, 
respectively, and reaches full current-sheet coverage at the maximum allowed refinement 
level 9 for t > 70.003. Current-density magnitude is color-shaded in the bottom panels, 
which are zoomed views of the initial magnetic null line (the red rectangle in the inset 
figure, top). 
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t = 0.0 s 


t = 147.64 s 




Figure 6. System configuration at the initial and final times t = 0.0 and 147.64 seconds. 
Field lines are shown in white. The background color scale shows the current strength 
along the perpendicular direction. The singular current system is seen clearly between the 
dislocated spine lines. 
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t = 0.0 s 


t = 100.1 s 

0.059 d 




t= 111.66 s 

0 .069 d 




t = 115.74 s 

0.072 d 



t = 147.64 s 

0.109 d 







Figure 7. Global and current-sheet structure at times t = {0.0, 100.1, 111.66, 115.74, 
147.64}. The fully developed current sheet has little curvature and retains an overall 
length L = 0.04 d. Current-density magnitude is color-shaded in each panel. 
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Figure 8. Time sequence showing the formation and ejection of a magnetic island (3D 
plasmoid cross-section), marked by the black arrows, in the perpendicular yz plane at x 
0. Plasmoid formation through linear growth and nonlinear saturation occurs between 
times t = 1 1 1.66 and 1 12.54, while plasmoid ejection occurs between times t = 1 13.08 
and 1 14.95. 
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Figure 9. Saturated shading shows the 3D structure and evolution of dense plasmoids 
(red) in magnetic islands along the entire strong current-density structure (purple ), in the 
xy plane viewed as from above. The white arrows mark the plasmoid tracked in cross- 
section in Figure 8. 
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t= 100.1 MMH^' 



Figure 10. Velocity magnitude color-shaded at the mid-plane at times t = {100.1, 1 1 1.66, 
1 15.74, 147.64} s. The reconnection jets are seen clearly as the purple shaded regions. 
Black contours are current-density magnitude, showing the approximate extent of the 
strong current sheet. White arrows mark the flow direction, showing the reconnection 
inflows from the side and the outflows at the ends of the current sheet. Note that the 
arrows are not scaled to the local velocity magnitude. 
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Figure 11. Evidence for a localized and transient Petschek-like geometry due to splitting 
of the current sheet at time t = 1 1 1.66 by plasmoid ejection. Top: Color shading of 
current density, red signifying strong currents (|J| > 5), with grid lines overlaid. Bottom: 
Color shading of plasma density, red signifying high densities (p > 2), with field lines 
overlaid. The hard bends in the field lines at the right-hand reconnection jet are Alfven 
waves associated with plasmoid ejection. 
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Figure 12. Large-scale periodic structure developing in the field lines along the third 
(initially axisymmetric) coordinate, exhibiting four wavelengths across the domain. 
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